home *** CD-ROM | disk | FTP | other *** search
/ Ham Radio 2000 #1 / Ham Radio 2000.iso / ham2000 / hf / dsp / dsp4tool / fft.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-02-19  |  3.3 KB  |  170 lines
  1. /*  FFT.C - Fast Fourier Transfom Pagage
  2.  *
  3.  *  Copyright (C) 1988-1994 by Alef Null. All rights reserved.
  4.  *  Author(s): J Vuori
  5.  *  Modification(s):
  6.  */
  7.  
  8. #include <math.h>
  9. #include "fft.h"
  10.  
  11.  
  12. /*
  13.  * Direct (Window)
  14.  *
  15.  * dest - pointer to complex array receiving the result
  16.  * src    - pointer to double array to be windowed
  17.  * n    - lenght of made transform (in form 2^m)
  18.  */
  19. void DirectWin(struct complex *dest, int *src, int n) {
  20.     struct complex *d;
  21.     int          *s;
  22.     int         i, m;
  23.  
  24.     m = (1 << n);
  25.     s = src; d = dest;
  26.     for(i = 0; i < m; i++) {
  27.     d->x = (double)*s;
  28.     d->y = 0.0;
  29.     s++; d++;
  30.     }
  31. }
  32.  
  33.  
  34. /*
  35.  * Hanning Window
  36.  *
  37.  * dest - pointer to complex array receiving the result
  38.  * src    - pointer to double array to be windowed
  39.  * n    - lenght of made transform (in form 2^m)
  40.  */
  41. void HanningWin(struct complex *dest, int *src, int n) {
  42.     struct complex *d;
  43.     int            *s;
  44.     int         i, m;
  45.  
  46.     m = (1 << n);
  47.     s = src; d = dest;
  48.     for(i = 0; i < m; i++) {
  49.     d->x = 0.5 * (1 - cos(2 * M_PI * (i + 0.5) / m)) * (double)*s;
  50.     d->y = 0.0;
  51.     s++; d++;
  52.     }
  53. }
  54.  
  55.  
  56. /*
  57.  * Hamming Window
  58.  *
  59.  * dest - pointer to complex array receiving the result
  60.  * src    - pointer to double array to be windowed
  61.  * n    - lenght of made transform (in form 2^m)
  62.  */
  63. void HammingWin(struct complex *dest, int *src, int n) {
  64.     struct complex *d;
  65.     int             *s;
  66.     int         i, m;
  67.  
  68.     m = (1 << n);
  69.     s = src; d = dest;
  70.     for(i = 0; i < m; i++) {
  71.     d->x = (0.54 - 0.46 * cos(2 * M_PI * (i + 0.5) / m)) * (double)*s;
  72.     d->y = 0.0;
  73.     s++; d++;
  74.     }
  75. }
  76.  
  77.  
  78. /*
  79.  * Blackman Window
  80.  *
  81.  * dest - pointer to complex array receiving the result
  82.  * src    - pointer to double array to be windowed
  83.  * n    - lenght of made transform (in form 2^m)
  84.  */
  85. void BlackmanWin(struct complex *dest, int *src, int n) {
  86.     struct complex *d;
  87.     int          *s;
  88.     int         i, m;
  89.  
  90.     m = (1 << n);
  91.     s = src; d = dest;
  92.     for(i = 0; i < m; i++) {
  93.     d->x = (0.42 - 0.5 * cos(2 * M_PI * (i + 0.5) / (m - 1))) * (double)*s;
  94.     d->y = 0.0;
  95.     s++; d++;
  96.     }
  97. }
  98.  
  99.  
  100. /*
  101.  * Fast Fourier Transform
  102.  *
  103.  * x - pointer to complex array, result replaces data that is here
  104.  * m - lenght of data (in form 2^m)
  105.  */
  106. void fft(struct complex x[], int m) {
  107.     register struct complex *a, *b;
  108.     double             w, dw,
  109.                  c, s,
  110.                  xt, yt;
  111.     int              i, j, k,
  112.                  n1, n2;
  113.  
  114.     n2 = (1 << m);
  115.  
  116.     for(k = 0; k < m; k++) {
  117.     /* handle all partitions */
  118.     n1 = n2;
  119.     n2 = n2 / 2;
  120.  
  121.     w  = 0.0;
  122.     dw = (2 * M_PI) / (double) n1;
  123.  
  124.     for(j = 0; j < n2; j++) {
  125.         /* handle all nodes */
  126.         c = cos(w);
  127.         s = sin(w);
  128.  
  129.         for(i = j; i < (1 << m); i += n1) {
  130.         /* basic butterfly (DIF) */
  131.         a = &x[i];
  132.         b = &x[i+n2];
  133.  
  134.         xt    = a->x - b->x;
  135.         a->x += b->x;
  136.         yt    = a->y - b->y;
  137.         a->y += b->y;
  138.  
  139.         b->x  = c * xt + s * yt;
  140.         b->y  = c * yt - s * xt;
  141.         }
  142.         w += dw;
  143.     }
  144.     }
  145. }
  146.  
  147.  
  148. /* Do bit-reversals on fft-result table
  149.  *
  150.  * i - direct index to result table
  151.  * n - lenght of made transform (in form 2^m)
  152.  */
  153. int Bit_Reverse(int i, int n) {
  154.     register int j, k;
  155.     int      m;
  156.  
  157.     m = (1 << n);
  158.     j = 0;
  159.  
  160.     for(k = 0; k < n; k++) {
  161.     m = m / 2;
  162.     if(i >= m) {
  163.         j += (1 << k);
  164.         i -= m;
  165.     }
  166.     }
  167.  
  168.     return(j);
  169. }
  170.